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Abstract 

We derive tight bounds on the cache misses for evaluation of explicit stencil oper- 
ators on structured grids. Our lower bound is based on the isoperimetrical property 
of the discrete octahedron. Our upper bound is based on a good surface to volume 
ratio of a parallelepiped spanned by a reduced basis of the interference lattice of a 
grid. Measurements show that our algorithm typically reduces the number of cache 
misses by a factor of three, relative to a compiler optimized code. We show that stencil 
calculations on grids whose interference lattice have a short vector feature abnormally 
high numbers of cache misses. We call such grids unfavorable and suggest to avoid 
these in computations by appropriate padding. By direct measurements on a MIPS 
R10000 processor we show a good correlation between abnormally high numbers of 
cache misses and unfavorable three-dimensional grids. 


1 Introduction 

On modern computers the gap between access times to cache and to global memory amounts 
to several orders of magnitude, and is growing. As a result, improvement in usage of the 
memory hierarchy has become a significant source of enhancing application performance. 
Well-organized data traffic may improve performance of a program, without changing the 
actual amount of computation, by reducing the time the processor stalls waiting for data. 
Both data location and access patterns affect the amount of data movement in the program, 
and the effectiveness of the cache. 

A number of techniques for improvements in usage of data caches have been developed 
in recent years. The techniques include improvements in data reuse (i.e. temporal locality) 
[3, 4, 5, 13], improvements in data locality (i.e. spatial locality) [13], and reductions in 
conflicts in data accesses [4. 5, 9, 10]. In practice, these techniques are implemented through 
code and data transformations such as array padding and loop unrolling, tiling, and fusing. 
Tight lower and upper bounds on memory hierarchy access complexity for FFT and matrix 
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multiplication algorithms arc given in [8]. However, questions concerning bounds on the 
number of cache misses and how closely current, optimization techniques approach those 
bounds for stencil operators remain open. 

In this paper we consider improvement of cache usgae through maximizing temporal 
locality in evaluations of explicit stencil operators on structured discretization grids. Our 
contribution if twofold. First, we prove lower and upper bounds on the number of cache 
misses for local operators on structured grids. Our lower bound (i.e. the number of un- 
avoidable misses) is based on the discrete isoperimetric theorem. Our upper bound (i.e. the 
achievable number of misses) is based on a cache fitting algorithm which utilizes a special 
basis of the grid interference lattice. As shown bv example, the lower bound can be achieved 
in some cases. The second contribution is the identification of grids unfavorable dimensions 
which cause significant increases in cache misses. We provide two characterizations of these 
unfavorable grids. The first one, derived experimentally, states that the product of all rele- 
vant grid dimensions is close to a multiple of half the cache size. The second characterization 
is that the grid interference lattice has a short vector. 

2 Cache model and definitions 

We consider a single-level, virtual-address-mapped, set-associative data cache memory, see 
[7]. The memory is organized in a sets of 2 lines of w words each. Hence, it can be 
characterized by the parameter triplet (a, z, w), and its size S equals a* z*w words. A cache 
with parameters (a, l,tn) is called fully associative, and with parameters (1, z.w) it is called 
direct-mapped. 

The cache memory is used as a temporary fast storage of words used for processing. A 
word at virtual address A is fetched into a {a(A), z(A),w(A)) cache location, where tu(A) = 
A mod w, z(A) = ( a/w ) mod z , and a(A) is determined according to a replacement policy 
(usually a variation of least recently used). The replacement policy is not important within 
the scope of this paper. 

If a word is fetched, then w — 1 neighboring words are fetched as well to fill the cache line 
completely. In practice, a, z, and w are often powers of 2 in order to simplify computation 
of the location in cache. For example, the MIPS R10000 processor for which we report some 
measurements in Section 6, has a cache with parameters (2,512,4), which makes S equal to 
4K double precision words, or 32KB. 

Our lower bound for the minimum number of cache misses that must be suffered during 
a stencil computation holds for any cache, including fully associative caches. The upper 
bound shows that a particular number of cache misses can be achieved by choosing a special 
sequence of computations. A cache miss is defined as a request for a word of data that is 
not present in the cache at the time of the request. A cache load is defined as an explicit 
request for a word of data for which no explicit request has been made previously (a cold 
load), or whose residence in the cache has expired because of a cache load of another word 
of data into the exact same location in the cache (a replacement load). The definitions of 
cold and replacement loads match those of cold and replacement cache misses, respectively 
[4], and if w equals 1 they completely coincide. 

If a piece of code features (p cache misses and p cache loads, it can easily be shown that 
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ft < wo. For a code with good spatial locality we typically have / 1 ~ wo. As can be shown by 
a simple example, no bound of the form 6 < r/t (r constant) can bo derived for arbitrary code 
segments, but if the code implements a noil-redundant stencil operation, we have o < | AT j // . 
where |A"| is the total number of points within the stencil. This is shown as follows. Let the 
stencil operation be written as q(x) = A’u(x), with x G fb Here id is the (not necessarily 
contiguous) point set on which array q is evaluated. Let 12 be the A'-extension of id. which 
is the point set on which u must be defined in order to compute q at all points of id. The 
total distinct number of elements of u used is |id|. The number of cache misses (j) does not 
exceed the total number of accesses to array u (may included repeated accesses to the same 
element), which equals |A||Q|, so [id j < |A'||12|. Consequently, w r e have the following interval 
inequality: |A’| -1 < ^ < w , which can be used to bound the number of cache misses in terms 
of the number of cache loads. 

3 A lower bound for cache loads for local operators 

In this section we consider the following problem: for a given d-dimensional structured grid 
and a local stencil operator A', how many cache loads have to be incurred in order to compute 
q = Ku, where q and u are two arrays defined on the grid. We will provide a lower bound p 
which asserts that, regardless of the order the grid points are visited for the computation of q , 
at least p cache loads have to occur. In the next section we provide a cache fitting algorithm 
for the computation of q whose number of cache loads closely approaches the lower bound. 

We use the following terminology to describe the operator K. The vectors ki,...,k 5 
defined such that q(x), the value of q at the grid point identified by the vector x, is a 
function of the values u(x + k t ),. . . ,u(x + k s ), are called stencil vectors. Locality of K means 
that the stencil vectors are contained in a cube {x| < r, i = 1 ,...,d} (r is called the 

radius of AT, and 2r + l its diameter). In this section we assume that K contains only the 
star stencil (i.e. the {0, ei, . . . , e<f, — ej, . . . , — e^} stencil). A lower bound for cache loads for 
the star stencil will give us a lower bound for any stencil containing it. 

Let q be computed in the AMnterior A of a rectangular region (a grid ) G. We assume 
that computation of q is performed in a pointwise fashion, that is, at any grid point the 
value of q is computed completely before computation of the value of q at another point is 
started. In order to compute the value of q at a grid point x, the values of u in neighbor 
points of x must be loaded into the cache (a point y is a neighbor of x if y — x is a stencil 
vector of A'). If x is a neighbor of y and u(y) has been loaded in cache to compute q( z) but 
is dropped from the cache before q(x) is computed, then u(y) must be reloaded, resulting in 
a replacement load associated with x. 

To estimate the number of elements, p, of array a that must be replaced, we choose a 
partition of A into a disjoint union of grid regions R t , with R = uf =1 A,, in such a wav that q 
is computed in all points of R t before it is computed at any point of A I+l , see Figure 1. Let 
B t j be the set of points in Rj which are neighbors of Aj. Since the star stencil is symmetrical, 
the B l} are neighbor points of Ajj. Because any point of B t] can have at most 2d neighbors 
in B }1 , we have the following inequalities: 

^\B lj \<\B ]l \<2d\B lj \. (1) 
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Figure 1: The boundaries Bi } of already computed values of q in a sequence of regions R l . 
Reloading of some values of u on the boundary of R$ results in at least max(|£? 3 i| + |R 32 | — S, 0) 
cache loads. 


For computation of q in R t we have to replace at least pi values of u , where p t equals 
max (e;;‘. \b,j\ — S, Oj . The total number of replaced values in the course of computing q 

on the entire grid will be at least p , where p equals B — kS, and B equals ^* =1 |Rq|. 

Summing all the terms \B l} \, taking into account Equation (1) and the fact that \Bu\ = 0, 
we get: 


B 


- k k 

4eei«, 


( 2 ) 


»=1 j=zl 


Let 6Ri be the exterior boundary of R u that is, all points neighbor to Ri not belonging 
to R t , and let 5 t be the subset of the grid boundary D having neighbors only in R t . Here, D 
is defined as G\R. Obviously, 5i O 5j = 0 if i ^ j, and ]T)j=i |jBq| > — |5 t | • 

Now t we choose the R t such that |(5Rj| = o > 8 dS, where a is specified below, and let u 


equal max|i?,|. Consequently, we have k > | R\jv = V/v, and 

k k 

p > - wo - kS = k (a - s ) ~ sE i*i a 7 s - i |Di 

1= 1 l— 1 

We subsequently choose a in such a way that 

a = \60(d,t)\ = j2 2k ( d k) {k-l)~ 8dS 


(3) 


(4) 


for some t, where 0(d, t) is the standard d-dimensional octahedron of radius t (see Appendix 
A). It follows from Equation 21, Appendix A, that t can be chosen in such a way that a is 
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loss than Sd(2d 4- 1)5. Now the value of u ran he estimated using the isoperimetric property 
of the octahedron (see again Appendix A), namely: // < \0{d,t)\. Hence, we find 


S n \SO(dJ)\ 

— > > ! 1 — — > ciS * ~ 

v ~ 8d(2d + l)v — Sd(2d + l)|0(rf, t)| — 


( 5 ) 


where c d equals \/{d{2d 4- l)2 d+2 ). This gives the following lower bound: 


P > 


V 


\50(d.t)\ 


1 


8d(2d+l) \0{d,t)\ 4 d 


\D\> Vc d S~ 


t^ d '• 


( 6 ) 


We also have V 4- \D\ = \G\ and \D\ < 2d\G\/l, where / is the smallest size of the grid. This 
gives the final lower bound fj. for the total number of elements of a to be read into the cache: 

P > V + p>V (l + c d S ~^ t - ^ > \G\ ^1 - — j — + (1 - y)c d 5“5 r >^ . (7) 

In general, assuming that the cache associativity a is larger than the diameter of the 
operator K, the order of this lower bound can not be improved, as shows the following 
example (remember that our lower bound is valid for a cache with any associativity, including 
a fully associative cache). Let the spatial extents of a two-dimensional grid be ri\ and ^ 2 , 
respectively, with riy equal to kS and ri 2 arbitrary, and perform calculations of the star stencil 
(i.e. r = 1) in the following order: 


do i = 0, k*a-l 
do j = 2, n 2 -l 

do il = max(2,l+i*(S/a)) , minCni-l, (i+l)*(S/a)) 
q(il» j) = u(il, j) + • • • 
end do 
end do 
end do 


Since rii equals kS, all values of q and u having the same value of the second index are 
mapped into the same cache location within a set. Since a exceeds 2r + 1, none of the values 
required for the computation of q will be replaced in the cache, except those at a distance 
r around the line defined by il = i*S/a. The total number of elements of u read into the 
cache for execution of this loop nest will therefore be nyn 2 + (n 2 — 2)2 r(ka — 1) — 4, which 
equals n l n- 2 ( 1 — 2/n { + 2a(l — 2/n 2 )/5). Similar examples in higher dimensions show that 
the order of our lower bound (Equation 7) can not be improved. 


4 An upper bound for cache loads for local operators. 
Cache fitting algorithm 

In order to obtain an upper bound we present a cache fitting algorithm which has a small 
number of replacements. We find a set P of cache conflict-free indices of u and calculate Ku 
at the points of P. Then we tile the index space of u with P to minimize the total number of 
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replacements. For the analysis wo assumo an eacho associativity of one. which is tho worst 
ease for replacement loads. 

Let L be a set in the index space of a having the same image in cache as the index 
(0, . . . , 0), Figure 2. L is a lattice in the sense that there is a generating set {b;}, i = 1, . . . , d, 
such that L is the set of grid points { (0, . . . , 0) 4- ^“ =1 x,b, | x, G Z }. We call this the 
interference lattice of a. It can be defined as the set of all vectors (q, • • ■ . i,i) such that 

i\ 4- n ih + n l n 2 i.i 4- • • • 4- n t • • • nd-Hd = 0 mod S . (8) 

In [4] this lattice is defined as the set of solutions to the cache miss equation. 

Let P be a fundamental parallelepiped of L*. For future reference we note that vol(P) = 
det L = S. The second equality follows form the fact that L has a basis {v ( } of the form: 

Vi^Pei, Vj = — mjei 4-e* ,2 < i < d,m l+l = • (9) 

j=i 

Obviously, the vectors satisfy Equation 8. Conversely, any vector satisfying Equation 8 can 
be represented as a linear combination of v lT . . . , v<*, with coefficients x * = 4 for k = 2, ... , d, 

and Xi = ( ii + m 2 i 2 4 F m d id)S~ x . x l is an integer number, since i\ + m 2 i ‘2 4 1- rn d id 

is divisible by S according to Equation 8. Since v^ . . . , v d are linearly independent vectors, 
they form a basis of the lattice. 

Let F be a face of P (see Figure 2), and let v be a basis vector of L such that P = 
{f 4- xv | f G F, 0 < x < 1}. Then shifts F 4- (k/g)v, k = . . . , -1,0, 1 , . . . contain all integer 
points of a pencil Q, with Q = {f 4- xv | f G F,x is any number} for an appropriate value 
of tjL The values of q at the points of Q can be computed without replacing reusable values 
of u except at a distance of r or less from the boundary of Q. Let hi, . . . ,h s be the signed 
projections along F of the stencil vectors of K onto v, and let h + and /i_ be the maximum 
and the minimum of the projections, respectively. We assume also that \h + — h_\/g < |v|a, 
meaning that the extent of P in the direction of v is big enough to allow to compute q on 
F without replacements. It may be impossible to satisfy this condition when the shortest 
vector in L is shorter than the diameter of K divided by the cache associativity. Lattices 
with short vectors are discussed in Section 6. The associated grids are called unfavorable 
grids. 

The Cache Fitting Algorithm for computing q is as follows (see Figure 2); here K(R) is 
the set of points where u must be available in order to compute q in all points of R (i.e. the 
A'-extension of R): 

set w = (1/g) v 

do Q = Qmin, Qmax 

determine face F inside pencil Q 

*A fundamental parallelepiped of a lattice L is a set of points x,b, 1 0 < x, < 1} for any basis {b, } 

of L. 

Wet / be a fundamental parallelepiped of the integer lattice in the subspace Y generated by F , and let 
e be an integer vector such that e and the basis of I generate Z d . Obviously, g must be chosen in such a 
way that vol((l/g)v. I) = vol (e, /) = 1. Hence, g = vol(v, I) = (l/|Fj)vol(v, F) = vol(P)/|F|, where |F| is 
the index of the lattice L fl Y in the integer lattice of Y. 
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do k = kmin, kmax 

load in cache all values of u inside /\'( F + k * w) 
compute q at F + k * w 
end do 
end do 

In this algorithm the parameters Qmin, Qmax, kmin and kmax are determined such that 
the scanning face F sweeps out the entire grid. Whenever a point is not contained in the 
grid, it is simply skipped in the nest. 

Since we defined the algorithm in such a way that the scanning face in the direction of v 
with step g passes through all integer points of Q. the values of q at all points inside Q will 
be computed. 



Figure 2: The Interference Lattice. Cache fitting set F + Arw, k E Z, sweeps across pencil 
Q in the direction of v. Only values of u at points at a distance r or less from the pencil 
boundaries /3i and /?2 will be replaced in the cache when K is evaluated inside of Q. 

Replacements misses can occur only at points at a distance of r or less from the boundaries 
of the pencils. For each of these points at most s replacement need to take place, where s, 
the size of the stencil, is defined by s = \K\ < (2 r + l) d . So the number of replacements will 
not exceed r(2r + l) d .4, where .4 is the total surface area of all pencils. 

To minimize .4 we choose P so that Q has a good surface to volume ratio. Let P be the 
fundamental parallelepiped of a reduced basis of L. A basis bi, . . . , b</ of a d-dimensional 
lattice L is called reduced if 

IJ||bi|| < c d det£r, ( 10 ) 

l 
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whore c d is a roust ant which depends only on r/L Let b t 1h' the shortest vector of the basis, 
and let the eccentricity e of f ile basis be defined by e = mnx( 1 1 b t M / 1 1 b 1 1 1 ) - If we define dP 
as tiie surface of P, we can derive an estimation for the surface-to-voiurne ratio of P: 


m 

det L 


< 


INI 


det L 


< '2c d ^ 


|bj| 


< c' 


llb.l! 




< ec' d S-* 


, 11 ) 


where we twice used the Hadamard inequality: It INI > detL, and the abovementioned 
fact that L = S. The constant c' d is defined by c' d = 2dc d . 

Since A does not exceed the surface area of all fundamental parallelepipeds covering 
the grid, the total number of these parallelepipeds (which equals \G\j detL) gives us: .4 < 
\dP\\G\j detT, so that the total number of replacements p can be bounded by: p < r(2r + 
l) d \dP\\G\/detL. This, combined with Equation 11, gives an upper bound for the total 
number of elements to be loaded into the cache in the cache fitting algorithm: 

p. < 1^1 + p < \G\ ^1 + ec" d S ^ , (12) 

where c d is defined by c d = r(2r + Vj d c' d . 

Note that if the shortest vector in the interference lattice has length (5/ f) l ^ d for some 
constant / it follows that e < fc d . To show this, we sort the basis vectors in Equation 10 in 

ascending order. Then it follows that j d ||bd|| < c d S, and hence e = < fc d . 

In Appendix B we show that there are grids w r hose interference lattices feature /’ s that 
are independent of S (provided that S is a prime power, which is true in most practical 
cases). For these lattices the relative gap between the upper bound (Equation 12) and the 
lower bound (Equation 7) of the previous section goes to zero as S increases. When the cache 
associativity exceeds the diameter of K, this gap can be closed. In that case a parallelepiped, 
built on a reduced basis of the interference lattice of the array indices with x d = 0, can be 
swept in the d th coordinate direction, similar to the example at the end of Section 3. In 
general, the cache fitting algorithm gives full cache utilization, in contrast to the algorithm 
for finding grid-aligned parallelepipeds devoid of interference lattice points, as proposed in 
[4]. See Table 2, [4], where the sizes of blocks without self interference are approximately 
20% smaller than S. 


5 Lower and upper bounds for multiple RHS arrays 

In this section we consider the case where there multiple arrays involved in the computation 
of q. Let p be the number of arrays (we call these the RHS arrays), all having the same sizes, 
and let the stencil of each RHS array include the star stencil. This means, in particular, 
that for each boundary point of any region R t (see Figure 1) values of p RHS arrays are 
necessary^ for computation of q in R x . Hence, we have to replace at least p t values, with 

* Every lattice has a reduced basis. There is a polynomial algorithm to find a reduced basis with a constant 
<71 = 2' <(d - 1 ^ 4 [11, Ch. 6.2]. 

§As in Section 3, we assume that computation of q is performed in a pointwise fashion. In this case 
elements of all RHS arrays have to be loaded into cache simultaneously, reducing the cache size effectively 
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p, = niiix{p(Y^j. , ! Bij | ) — 5.0) values of RHS arrays. Now we can repeat the arguments of 
Section 3. with 1 1 ’ | and |G'| replaced bv p|l’| and p\G\. respectively, and 5 replaced by \S/p\. 
to obtain the following lower bound fur the number of cache loads for stencil computations 
with p RHS arrays: 


p > p\V\ + p > p\V\ 1 + Q 


> p\G\ 1 - 


2d- 1 

1 


i 

~ 21 

> 

+ 1 — r I c <i 


'S' 

P 


1 

' d- 1 


(13) 


In order to obtain an upper bound for cache loads for calculations with p RHS arrays, 
we assume that we are free to choose relative array offsets. Our upper bound is valid on 
the assumption that the stencil diameter divided by the cache associativity is smaller than 
the length of the longest lattice basis vector divided by p. Consider a stripwise tiling of 
the fundamental parallelepiped P for the lattice L, see Figure 3. Each tile P t has the same 
size and shape. The size is determined by considering the longest edge vector v in the 
fundamental parallelepiped and dividing it into p equal pieces of size [(5/|F|)/p]||v||. so that 
each tile contains |F|[(5/|F|)/p] integer points, where |F| is the number of integer points in 
the face. The remainder part of the tiling is indicated by the shaded area. The reason why 
the longest edge vector is selected for subdivision is as follows. Since w r e use a reduced basis, 
the smallest angle between v and F is bounded from below, so the parallelepiped is always 
close to orthogonal. Therefore, subdividing the longest edge leads to tiles with the largest 
inscribed sphere, and thus the largest difference stencil fitting inside the tile. 

Let {Pi} be the parallelepipeds of the tiling, and let s* be the address offset of F, relative 
to Pi (corresponding to the same RHS array). We assign one parallelepiped to each RHS 
array and choose starting addresses of the arrays, addr*, in such a way that images of tiles 
Pi in the cache do not overlap: addr* = addr t + rriiS + s,, where m\ = St = 0, and 
m l = m,_i + i = 2 ,p. Sweeping through the pencil by units of tile P\ in the 

direction of v we can compute Ku without any cache conflicts, except on the boundary of 
the pencils. The number of replacement loads of this algorithm can be estimated similarly 
to the number of replacement loads of a single-array algorithm, taking into account that for 
calculation of a value u at any point values of all p RHS arrays in the neighbor points may 
have to be in cache, thus reducing the effective cache size to [5/pj: 


(14) 

where c " d is a constant which depends only on d, and e is the eccentricity of L. 

by a factor of p. Non-pointwise computations may be performed if the operator K is separable, in the sense 
that it can be written as A'(ui, u-z , . . . , u p ) = K\ (ui , . . . K p (u p ) . . . ). In the case of separability of A 

the stencil operation can be split into a succession of independent operations, each involving an intermediate 
value of q and one RHS array. This would not require to load all p RHS arrays in cache at each point. 
Instead, it would suffice to write intermediate values of q into main memory, and then load them back into 
cache for completion of the computations. This results in a larger effective cache size, but more data to be 
loaded, so splitting the operation need not improve the total number of loads. 


P < p\G\ + P < p\G\ fl +ecj 


5 

IP 


9 




Figure 3: Tiling of a fundamental parallelepiped of a reduced basis of the lattice L. We 
assume that \h + — 1 < a|v/p| ( a is the cache associativity). The tiling effectively reduces 

the size of the parallelepiped by a factor of at most 2 p (since x/p > \_xjp\ > x/(2p)), and 
increases the cost of a replacement in the cache per point of the boundary of the pencil by 
at most a factor of p, since elements of all p RHS arrays will be replaced at the same time. 


6 Unfavorable array sizes 

We have implemented our cache fitting algorithm and compared its actually measured num- 
ber of cache misses with those of the compiler-optimized code for the corresponding naturally 
ordered loop nest on a MIPS R10000 processor (SGI Origin 2000). For comparison we chose 
a second order difference operator (the common 13-point star stencil) an a test set includ- 
ing three-dimensional grids of sizes 40 < ri\ < 100, n ,2 = 91, and n 3 = 100 (the value of 
the second dimension was chosen to show a typical picture; that of the third dimension is 
irrelevant). A plot of measured cache misses for both codes is shown in Figure 4. The pro- 
gram was compiled with options “-03 -LN0 :prefetch=0,” using the MIPSpro f77 compiler, 
version 7.3.1.1m. The prefetch flag disables the prefetching compiler optimization. Without 
this option the number of cache misses increases significantly, because the compiler does 
aggressive prefetching to try to reduce execution time. 

The upper bounds for the cache misses from the previous sections would suggest that the 
number of replacement cache misses will increase in the cases where the interference lattice 
has a very short vector. Very short means that the length is smaller than the diameter of the 
operator divided by the cache associativity. In this case the self interference would increase 
significantly. This result suggests how to pad arrays to improve cache performance: the 
padding should be organized in such a way that the shortest vector in the lattice is not too 
short, though short enough to minimize the number of pencils (large index of scanning face 
F). The sweeping is organized such that pencils are as wide as possible (i.e. the smallest 
total number of pencils), while avoiding — in the case of multiple RHS arrays — tiles that are 
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second array dimension (n^ 


Figure 4: Plot of measured cache misses for 40 < n\ < 100, n 2 = 91 for 13-point star stencil. 
The top line corresponds to the naturally ordered nest, optimized by the SGI Fortran com- 
piler. The bottom line corresponds to our cache fitting algorithm. A typical ratio between 
the two is 3.5. The large fluctuations correspond to grids with short lattice vectors (n t = 45 
and = 90 yield shortest vectors (1,0,1) and (2,0,1), respectively). The fluctuations of 
cache misses of the cache fitting algorithm for such grids can be so big that their cache misses 
become more numerous than for the compiler-optimized nest. 


thinner than the diameter of the stencil operator divided by the cache associativity. 

To demonstrate these unfavorable grids we again choose the second order stencil and 
force computations in the nest to follow the natural orderV Figure 5a shows the correlation 
between spikes in the number of cache misses and the presence of a very short vector in the 
lattice. We call these lattices unfavorable for cache utilization. Arrays having such lattices 
should be avoided on the target machine. When the shortest vector of the interference lattice 
is shorter than the diameter of the operator, the number of cache misses sharply increases. 
The application developer should avoid such unfavorable array sizes, and compilers should 
avoid the sizes using appropriate padding of array dimensions. Note that similar unfavorable 
cache effects have been mentioned in [1], 

^This forcing is accomplished by introducing a dependence through a Fortran subroutine that performs 
a circular shift of its arguments 
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A) «/ — B) 




Figure 5: Plot A shows measured fluctuations of cache misses (above 15% of the upper 
bound). Plot B shows the interference lattices with short (less than 8 in the L\ norm) 
vectors. Array sizes are 40 < ni,n 2 < 100. The plots can be fitted well by hyperbolae 
defined by nin 2 = | kS,k = 1,2, 3, 4, meaning that arrays with unfavorable size are those 
whose 2 -slices are (close to) multiples of half the cache size. The horizontal line in Plot A 
shows the position of the graph from Figure 4. 


7 Conclusions and future work 

We have demonstrated tight lower and upper bounds for cache misses for calculations of an 
explicit operator K on a structured grid. Our lower bound is valid in the general case of 
fully associative caches, and is based on a discrete isoperimetric theorem. Our upper bound 
is based on a cache fitting algorithm which uses the fundamental parallelepiped of a special 
basis of the interference lattice to fit the data in the cache. The upper bound assumes that 
the shortest vector in the interference lattice is not too short. We have shown that there are 
grids whose interference lattices have this property. We have also shown that the presence 
of a very short vector in the lattice correlates with fluctuations of actual cache misses for 
calculation of a second order explicit operator on three-dimensional grids. The fluctuations 
occur on grids with unfavorable sizes, i.e. on those whose product of the first two dimensions 
is (close to) a multiple of half the cache size. 

Our results can be extended straightforwardly to implicit stencil computations (i.e. those 
of the form q <— K (q)) when the problem has a one-dimensional data dependence. Such a 
data dependence exists if computations of q at grid points can take place in an arbitrary 
order, except that there is a single index i for which q(x \, . . . , i , . . . , x d) must be evaluated 
before q(x i , . . . , i 4- «, . . . , x^) can be calculated (the constant c* is either +1 or -1). Clearly, 
the lower bound is not affected by the implicitness of K. The previously derived upper bound 
can still be achieved by prescribing the proper visit order of points within each parallelepiped, 
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of the scanning face direction within each pencil (positive or negative sweep direction), and 
of the visit order of subsequent pencils. This is always possible for a one-dimensional data 
dependency. 

Our results can also be extended to arrays that store more than one word per grid point 
(tensor arrays). The lower bound of Section 3 for operations with multiple right hand sides 
immediately applies to tensor arrays. The upper bound of that section also applies, provided 
the tensor components can be stored as independent subarrays. 

In a future study we plan to extend the results of this paper to more general implicit 
operators, to operators on unstructured grids, and to tensor arrays with restricted storage 
models. We intend to study more closely the dependence of cache misses on the size of the 
operator's stencil. We also plan to enhance the presented results by taking into account 
a secondary cache and TLB, and to formulate bounds for cache misses more directly than 
through the determination of cache loads. 


Appendix A: The simplex and the octahedron 


In this section vve list some basic facts on the number of integer points in the octahedron 
and simplex. The standard octahedron is defined as: 


0(d, t) = l x G Z d | |xj| < t 


(15) 


i — l 


and the standard simplex as: 

S(d, t) = /x G Z d | 0 < Xu • • • >Xd , ^ |xj| < t 


(16) 


2=1 


If we consider sections of the octahedron by planes X\ = k, k = —t, . . . , t, then for the 
number of integer points in the octahedron we get the following recurrence relation: 


e-i 

1 0(d, i)| = 1 0(d - 1, f)| + 2 £ 1 0(d - 1, k)\ . (17) 

i t=o 

This relation can be used to prove that 

|0(<U)l = E 2 ‘QG) <18) 

fc=0 \ / \ / 

and that 

I SO(d, t - 1)1 = \0(d, t ) - 0(d, ( - Dl = E 2 ‘ (t) (1 1 1) ' (19) 

Also, the relation 

|<iO(d,t)| = \SO{d,t - l)\ + \50{d - l,t)\ + \SO{d - l,t - 1)| (20) 
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shows that 

\<)(.){,L 01 < (2d + 1)| 6()(dJ - 1)| . (21) 

For the number of integer points in the simplex we have the following recurrence relation: 
\S(<L 0| = \S(d — 1, 01 + \S(d,t - 1)| . (22) 


This can be used to prove that cf. [6], Table 169, see also [12], Section 5: 


iswoi-E 


A:— 1 



(23) 


From Equations 18 and 21 it follows that \0(djt)\ < 2 d |5(d,0l- Also, since SO(d,t — 1) 
contains at least two nonoverlapping simplices S(d — 1,0 and can be covered by 2 d such 
simplices, we see that 

2|S(rf - 1, 01 < | SO(d, t - 1)| < 2 d \ S(d -1,01, d> 2. (24) 


Hence, if | S(d — 1,01 equals S, we have for d > 2: 

jjgtofll > > giM! = 2 -*fi /\ * y 1 > 2~ d+i 

\0{d,t)\ ~ \0(d, 0| - 2 d \S(d, 0| V d) ~ 


(25) 


since from Equation 23 it follows that if |S(d — 1,01 does not exceed S, then 1 + t/d does 
not exceed S l ^ d ~ l K 

The isoperimetric inequality [12], Theorem 2, asserts that the size of the boundary of a 
subset R in Z d is at least as big as the size of the standard sphere that contains |i?| points". 
It is easy to see that any standard sphere is sandwiched between two standard octahedrons 
whose radii differ by one. Since the octahedron has the largest volume for a given fixed-size 
boundary, Inequality 25 is true for any lattice body with a boundary of size S. 


Appendix B: The existence of grids with favorable lat- 
tices 

In order to prove that for every cache of size S = p n , where p is a prime number, there are 
grids with interference lattices whose shortest vector has a length l greater than (5/ Z) 1 ^, 
with / independent of 5, we show: 

a. For every dimensionality d there exists a lattice L of the same dimension whose basis 
has the form given in Equation 9 (Section 4), and whose shortest vector is sufficiently 
long, and 

b. a grid can be constructed that has L as its interference lattice. 

^The standard sphere, defined in [12], is the integer point set of minimal surface area for any given number 
of interior points. 
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Corollary: Since grids with dimensions n t + k t S. i — <1 have the same interference 

lattice for any non-negative integers A-,, any grid can he embedded in a favorable larger grid. 


Proof: 

Step a: Let a lattice L have a basis of the form of Equation 9. Any lattice vector x, which 
includes all basis vectors of L, with norm at most l must be a solution of the 
following system of inequalities: 

|x t | < / , i = 2, . . . , d 
|5xi + 7712 X 2 + • • • + m d x d | < l 



Existence of a solution to this system is equivalent to that of the system 


N 

< 1 . 

,7=2,... 

,d 

) 

m 2 

+ ■ • 

m d 

< 

1 \ 

ly l2 


S > 


where ||^|| is the distance from 2 to the nearest integer number. Theorem VIII, Ch. 1 
[2] states that there are real numbers /i 2 , . . . , p d , and a constant c^" depending only on 
d, such that 


\\H 2 X 2 + • • • + fi d Xd\\ > 



(28) 


for all nonzero x satisfying |x;| < l , i = 2, . . . , d. 

If we choose the nonzero integers m x in such a way that |mj — S/q| < 2 for i = 2, . . . , d, 
then we get 


m 2 m d 

-JX 2 + • ■ • + -jx d 


> 



(d-C s 




(29) 


which shows that Equation 26 has no integer solutions if / < (c'^/(Sd)) l ^ d . Hence, / in 
Section 4 can be chosen as: / = d/c' d ", and the lattice with the basis given by Equation 
9 has a reduced basis with eccentricity depending only on d. 

Step b: In order to find a grid whose interference lattice is L, we first sort the m, in order of 
increasing gcd(m,,5). Since we assume that S = p n , where p is prime, we know that 
gcd(m,,5) divides gcd(m [+ i , S). and the appropriate grid dimensions n, can be found 
directly by solving the congruencies (n l m l — m 1+1 ) mod 5 = 0. 
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